#Librerias a usar
import numpy as np # Libreria numerica
import matplotlib.pyplot as plt #Libreria matemática
from scipy.io.wavfile import read # libreria para lectura de archivos de audio
from scipy.io.wavfile import write #libreria para crear el audio
from scipy.io import wavfile #Libreria para cargar el audio
from IPython.display import Audio # para escuchar la senal
import scipy.signal as sp #Libreria para señales
import sympy as sy
Frecuencia de corte 1: 3700 Hz. Frecuencia de corte 2: 4300 Hz. N=13. Fs=13000 Hz
fc1_f1=3700 # Frecuencia de corte 1
fc2_f1=4300 # Frecuencia de corte 2
N_f1=13 # orden del filtro
Fs_f1=13000 # Frecuencia de muestreo
#Se calcula M par poder saber cada cuanto se toma una muestra.
M_f1 = int(Fs_f1/N_f1) # Es para los pasos de k-> k(1khz)
K_f1 =1 # limite de la sumatoria
# Se realiza el siguiente condicional para cuando N sea par o impar, dado que la formula cambia segun corresponda.
if N_f1%2 ==0:
K_f1=int(N_f1/2 -1)# para un N par
else:
K_f1=int((N_f1-1)/2) #Para un N impar
# Se calcula el factor de simetría,el cual nos sirve para definir como se distribuyen las componentes con relación a su punto medio
alfa_f1 = int((N_f1-1)/2) #Constante
hk_f1 = [] #Coeficientes de mi ventana h para encontrar la ecuacion
n_f1 =np.arange(0,K_f1+1)
# Lleno el vector de coeficientes hk para multiplicar
for i in range(K_f1+1):
if i*M_f1 < fc1_f1 or i*M_f1>fc2_f1:
hk_f1.append(0)
else:
hk_f1.append(1)
#se crea una funcion symbolica
h_f1 = 1
n_f1 = sy.Symbol("n")
for i in range(1,int(K_f1)):
h_f1+=2*hk_f1[i]*sy.cos((2*np.pi*i*(n_f1-alfa_f1))/N_f1)
h_f1=(1/N_f1)*h_f1
#Array para sustituir valores de n de 0 a N/2
hn_f1 = np.array([])
for i in range(N_f1):
#Reemplazo la n completa para obtener los coeficientes del filtro
hn_f1= np.append(hn_f1,float(h_f1.subs(n_f1,i)))
w_f1,Hh_f1 = sp.freqz(hn_f1,1, whole = True, worN=1024) #Respuesta en frecuencia
for i in range(N_f1//2):
print("h({})".format(i),"= h({})".format(N_f1-i-1)," = ",np.round(hn_f1[i],5))
print("h(6) =",np.round(hn_f1[6],5))
h(0) = h(12) = 0.16432 h(1) = h(11) = -0.07245 h(2) = h(10) = 0.09547 h(3) = h(9) = 0.21315 h(4) = h(8) = -0.03823 h(5) = h(7) = 0.02237 h(6) = 0.23077
#Dominio de tiempo discreto
plt.figure(figsize=(15,7))
plt.stem(hn_f1)
plt.title("Respuesta al impulso del filtro",fontsize=20)
plt.xlabel("Muestras",fontsize=15)
plt.ylabel("Amplitud",fontsize=15)
plt.grid()
plt.show()
#Dominio de la frecuencia
plt.figure(figsize=(15,7))
plt.title('Respuesta en frecuencia del filtro pasa-banda por Muestreo en frecuencia', fontsize=20)
plt.plot((w_f1-np.pi)*Fs_f1/(2*np.pi), np.abs(np.fft.fftshift(Hh_f1)), label='$Filtro$')
plt.vlines(fc1_f1, 0, 1.2, color='r', lw=2, linestyle='--',label='$F_{c1}$')
plt.vlines(fc2_f1, 0, 1.2, color='g', lw=2, linestyle='--',label='$F_{c2}$')
plt.vlines(3000, 0, 1, color='b', lw=2, linestyle='--',label='$K=3$')
plt.vlines(5000, 0, 1, color='b', lw=2, linestyle='--',label='$K=5$')
plt.xlabel("$Frecuencia [Hz]$", fontsize=15)
plt.ylabel("$|H2(\omega)|$", fontsize=15)
plt.legend()
plt.grid()
plt.show()
def Audio_tiempo(t,y,audio,filtrado):
"""
Funcion para graficar el audio en el dominio del tiempo, tiene como entrada:
- t: Vector del eje x (tiempo)
- y: Vector para graficar los valores del eje y ()
- audio: el número del audio que estoy procesando para poder llevar el control de manera organizada
- filtrado: una variable booleana para verificar si se filtro o no el audio e imprimir el titulo
"""
plt.figure(figsize = (15,4))
if filtrado:
plt.plot(t,y,color="#746AB0")
plt.title("Señal filtrada del Audio {}".format(audio), fontsize=20)
else:
plt.plot(t,y,color="#E83845")
plt.title("Señal original del Audio {}".format(audio), fontsize=20)
plt.xlabel("Tiempo (s)", fontsize=15)
plt.ylabel("Amplitud", fontsize=15)
plt.grid()
plt.show() # Mostrar la grafica
def Audio_frecuencia(freq,xfreq,audio):
"""
Funcion para graficar el espectro de un audio original, tiene como entrada
- freq: vector de frecuencias
- xfrec: vector del espectro
- audio: número de audio para impresión en el titulo
"""
plt.figure(figsize=(15,7)) #Tamaño de la figura
plt.plot(freq, np.abs(xfreq),c="#ff7700")
plt.xlabel("Hz",fontsize=15)
plt.ylabel("$|H(\omega)|$",fontsize=15)
plt.title("Espectro original del Audio {}".format(audio), fontsize=20)
plt.xlim([-6000,6000])
plt.grid()
plt.show()
def G_Audio_Filtrado_freq(freq_fm,xfreq_fm,audio,fc1,fc2,):
"""
Funcion para graficar el espectro de un audio con filtro pasa-banda
con muestreo en frecuencia, tiene como entrada:
- freq_fm: vector de frecuencias filtradas
- xfrec_fm: vector del espectro filtrado
- fc1: Frecuencia de corte inferior
- Fc2: frecuencia de corte supereior
- audio: número de audio para impresión en el titulo
"""
plt.figure(figsize=(15,7)) #Tamaño de la figura
plt.plot(freq_fm, np.abs(xfreq_fm)/max(abs(xfreq_fm)),c="#388E8E")
plt.vlines(fc1, 0, 1, color='r', lw=2, linestyle='--',label='$F_{c1}$')
plt.vlines(fc2, 0, 1, color='g', lw=2, linestyle='--',label='$F_{c2}$')
plt.vlines(3000, 0, 1, color='b', lw=2, linestyle='--',label='$K=3$')
plt.vlines(5000, 0, 1, color='b', lw=2, linestyle='--',label='$K=5$')
plt.xlabel("Hz",fontsize=15)
plt.ylabel("$|H(\omega)|$",fontsize=15)
plt.title("Espectro filtrado del Audio {}".format(audio), fontsize=20)
plt.xlim([-6000,6000])
plt.legend()
plt.grid()
plt.show()
# Primero encontramos el ripple para determinar que tipo de enventanado se utiliza
rp_f2 =0.004 #ripple dado
rpdb_f2=20*np.log10(rp_f2) # ripple en db
print("El ripple en db:",rpdb_f2)
El ripple en db: -47.95880017344075
Y como es un filtro pasa bajas, entonces:
#Diseño del filtro y encontrar orden
fm_f2 =200 # Ancho de banda de transición
Fc_f2=600 # Frecuencia de corte
Fs_f2=12000 #Frecuencia de muestreo
M_f2 = (6*Fs_f2)/fm_f2 #Orden del filtro
wc_f2 = (2*np.pi*Fc_f2)/Fs_f2 #Frecuencia de corte normalizada en radianes
if(M_f2%2 == 0): #Lo hace impar
M_f2 = M_f2+1
M_f2= int(M_f2)
print('Orden del filtro Pasa-bajos: ', M_f2)
Orden del filtro Pasa-bajos: 361
#Creación de a ventana
n_f2 = np.arange(-(M_f2-1)/2,M_f2/2) #Vector de muestras
hi_f2 = (wc_f2/np.pi)*np.sinc((wc_f2*n_f2)/np.pi) #Respuesta del filtro ideal
win_f2 = 0.54-0.46*np.cos(2*np.pi*np.arange(len(n_f2))/(M_f2-1)) #ventana hamming
h_f2 = hi_f2*win_f2 #Multiplico la respuesta ideal por la ventana
#print(h)
w_sv_f2,Hh_sv_f2 = sp.freqz(hi_f2,1, whole=True, worN=2048) #Respuesta en frecuencia del filtro sin eventanado
w_f2, Hh_f2 = sp.freqz(h_f2,1, whole=True, worN=2048) #Respuesta en frecuencia del filtro enventanado
print("los 5 primeros y ultimos")
for i in range(0,6):
print("win({})".format(i),"= win({})".format(M_f2-i-1)," = ",win_f2[i])
print("h({})".format(i),"= h({})".format(M_f2-i-1)," = ",h_f2[i])
los 5 primeros y ultimos win(0) = win(360) = 0.08000000000000002 h(0) = h(360) = -3.1185374660155015e-19 win(1) = win(359) = 0.08007006022806001 h(1) = h(359) = -4.399968989045795e-05 win(2) = win(358) = 0.08028021957121595 h(2) = h(358) = -8.438346641468691e-05 win(3) = win(357) = 0.08063041401289606 h(3) = h(357) = -0.00011730955714971217 win(4) = win(356) = 0.08112053688048088 h(4) = h(356) = -0.0001395322512564273 win(5) = win(355) = 0.08175043887779704 h(5) = h(355) = -0.00014869698796952295
n_f2=np.arange(0,M_f2,1) # tamaño vector temporal
plt.figure(figsize=(15,7))
plt.stem(n_f2,hi_f2)
plt.title("Respuesta del filtro ideal sin enventanado",fontsize=20)
plt.xlabel("Muestras",fontsize=15)
plt.ylabel("Amplitud",fontsize=15)
plt.grid()
plt.show()
plt.figure(figsize=(15,7))
plt.plot(n_f2,win_f2)
plt.title("Ventana Hamming",fontsize=20)
plt.xlabel("Muestras",fontsize=15)
plt.ylabel("Amplitud",fontsize=15)
plt.grid()
plt.show()
plt.figure(figsize=(15,7))
plt.stem(n_f2,h_f2)
plt.title("Respuesta al impulso con enventanado",fontsize=20)
plt.xlabel("Muestras",fontsize=15)
plt.ylabel("Amplitud",fontsize=15)
plt.grid()
plt.show()
plt.figure(figsize=(15,7))
plt.title('Respuesta en frecuencia del Pasa-bajos sin enventanado', fontsize=18)
plt.plot((w_sv_f2 -np.pi)*Fs_f2/(2*np.pi), np.abs(np.fft.fftshift(Hh_sv_f2)), label='$Filtro$')
plt.vlines(Fc_f2, 0, 1, color='r', lw=2, linestyle='--',label='$F_c$')
plt.vlines(Fc_f2+fm_f2/2, 0, 1, color='g', lw=2, linestyle='--',label='$F_c+fm/2$')
plt.vlines(Fc_f2-fm_f2/2, 0, 1, color='g', lw=2, linestyle='--',label='$F_c-fm/2$')
plt.xlim([-0,2500])
plt.xlabel("$Frecuencia [Hz]$", fontsize=18)
plt.ylabel("$|H2(\omega)|$", fontsize=18)
plt.legend()
plt.grid()
plt.show()
plt.figure(figsize=(15,7))
plt.title('Respuesta en frecuencia del Pasa-bajos con enventanado', fontsize=18)
plt.plot((w_f2-np.pi)*Fs_f2/(2*np.pi), np.abs(np.fft.fftshift(Hh_f2)), label='$Filtro$')
plt.vlines(Fc_f2, 0, 1, color='r', lw=2, linestyle='--',label='$F_c$')
plt.vlines(Fc_f2+fm_f2/2, 0, 1, color='g', lw=2, linestyle='--',label='$F_c+fm/2$')
plt.vlines(Fc_f2-fm_f2/2, 0, 1, color='g', lw=2, linestyle='--',label='$F_c-fm/2$')
plt.xlim([-0,2500])
plt.xlabel("$Frecuencia [Hz]$", fontsize=18)
plt.ylabel("$|H2(\omega)|$", fontsize=18)
plt.legend()
plt.grid()
plt.show()
def filtro_enventanado(freq,xfreq,audio,fc,fm):
"""
Funcion para graficar el espectro de un audio con filtro pasa-baja
con muestreo en frecuencia, tiene como entrada:
- freq: vector de frecuencias filtradas
- xfre: vector del espectro filtrado
- fc: Frecuencia de corte del filtro
- audio: número de audio para impresión en el titulo
"""
plt.figure(figsize=(15,7)) #Tamaño de la figura
plt.plot(freq, np.abs(xfreq)/max(abs(xfreq)),c="#388E8E")
plt.vlines(fc, 0, 1, color='r', lw=2, linestyle='--', label='$F_c$')
plt.vlines(fc+fm/2, 0, 1, color='m', lw=2, linestyle='--', label='$F_c+f_m/2$')
plt.vlines(fc-fm/2, 0, 1, color='m', lw=2, linestyle='--', label='$F_c-f_m/2$')
plt.xlabel("Hz",fontsize=15)
plt.ylabel("$|H(\omega)|$",fontsize=15)
plt.title("Espectro filtrado del Audio {}".format(audio), fontsize=20)
plt.xlim([-2000,2000])
plt.legend()
plt.grid()
plt.show()
# Señal original de Audio
# cargar el archivo de audio estéreo
rate_f1, data_f1 = wavfile.read('naruto.wav')
left_channel_f1 = data_f1[:, 0]
norm1_f1 =left_channel_f1/np.max(np.abs(left_channel_f1))
t1_f1=np.arange(0, len(data_f1)/rate_f1, 1.0/rate_f1) # Vector de tiempo
print(rate_f1)
44100
Audio_tiempo(t1_f1,norm1_f1,1,False)
Audio(norm1_f1, rate=rate_f1)
# Espectro del audio 1 original
xfreq1_f1 = np.fft.fft(left_channel_f1) #Transformada Rapida de Fourier del audio
n_1_f1 = np.size(xfreq1_f1) #Vector de Frecuencias
freq1_f1 = np.arange(-n_1_f1/2, n_1_f1/2)*rate_f1/n_1_f1#Vector de Frecuencias = np.size(xfreq1)
xfreq1_f1 = np.fft.fftshift(xfreq1_f1) #organizar el vector del espectro:
Audio_frecuencia(freq1_f1,np.abs(xfreq1_f1)/max(abs(xfreq1_f1)),1)
# Aplicando primer filtro a señal de audio
y_1_f1 = np.convolve(hn_f1, norm1_f1, mode='same') #Filtra la señal requerida
Y_1_f1 = np.fft.fft(y_1_f1) #FFT de la señal filtrada para mostrarla más comodamente
freq_Y_f1 = np.fft.fftfreq(Y_1_f1.size)*Fs_f1 #Vector de Frecuencias
#Se grafica el audio de la señal filtrada en nuestro primer filtro
filtrado=True
Audio_tiempo(t1_f1,y_1_f1,1,True)
Audio(y_1_f1, rate=rate_f1) # Reproductor señal de audio filtrado
#Transformada Rapida de Fourier del audio
xfreq1_f1 = np.fft.fft(y_1_f1)
#Vector de Frecuencias
n_1_f1 = np.size(xfreq1_f1)
#Vector de Frecuencias = np.size(xfreq1)
freq1_f1= np.arange(-n_1_f1/2, n_1_f1/2)*rate_f1/n_1_f1
#organizar el vector del espectro:
xfreq1_f1 = np.fft.fftshift(xfreq1_f1)
#Normalizar amplitud del vector
# Grafica
G_Audio_Filtrado_freq(freq1_f1,xfreq1_f1,1,fc1_f1,fc2_f1,)
# Aplicando segundo filtro a señal de audio
y_1_f2 = np.convolve(h_f2, y_1_f1, mode='same') #Filtra la señal requerida
Y_1_f2 = np.fft.fft(y_1_f2) #FFT de la señal filtrada para mostrarla más comodamente
freq_Y_f2 = np.fft.fftfreq(Y_1_f2.size)*Fs_f2 #Vector de Frecuencias
#Grafica
#Se grafica el audio en el tiempo ya filtrado
filtrado=True
Audio_tiempo(t1_f1,y_1_f2,2,filtrado)
Audio(y_1_f2, rate=rate_f1) # Reproductor señal de audio filtrado
#Transformada Rapida de Fourier del audio
xfreq1_f2 = np.fft.fft(y_1_f2)
#Vector de Frecuencias
n_1_f2 = np.size(xfreq1_f2)
#Vector de Frecuencias = np.size(xfreq1)
freq1_f2 = np.arange(-n_1_f2/2, n_1_f2/2)*rate_f1/n_1_f2
#organizar el vector del espectro:
xfreq1_f2 = np.fft.fftshift(xfreq1_f2)
filtro_enventanado(freq1_f2,xfreq1_f2,2,Fc_f2,fm_f2)